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Large birds and glider pilots commonly use updrafts caused by convection in the lower 
atmosphere to extend flight duration, increase cross-country speed, improve range, or 
simply to conserve energy. Uninhabited air vehicles may also have the ability to exploit 
updrafts to improve performance. An updraft model was developed at NASA Dryden Flight 
Research Center (Edwards, California) to investigate the use of convective lift for 
uninhabited air vehicles in desert regions. Balloon and surface measurements obtained at the 
National Oceanic and Atmospheric Administration Surface Radiation station (Desert Rock, 
Nevada) enabled the model development. The data were used to create a statistical 
representation of the convective velocity scale, w , and the convective mixing-layer thickness, 
Zi- These parameters were then used to determine updraft size, vertical velocity profile, 
spacing, and maximum height. This paper gives a complete description of the updraft model 
and its derivation. Computer code for running the model is also given in conjunction with a 
check case for model verification. 
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= specific heat of dry air, J/kg°K 
= saturated vapor pressure, mb 
= vapor pressure, mb 
= gravitational constant, m/s 2 
= height above ground, m 
= ground measurement sample index 
= shape constant, unitless 
= straight line length, m 
= maximum 

= ground measurement sample index at sunset 
= number of updrafts inside a given area 
= number of updrafts encountered along straight line 
= surface pressure, mb 
= reference pressure, mb 
= probability density function 
= heat into ground, W/m 2 
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sensible heat flux, W/m 2 

kinematic sensible heat flux, °K*m/s 

surface virtual potential temperature flux, W/m 2 

net radiation at surface, W/m 2 

distance to updraft center, m 

updraft core radius, m 

updraft outer radius, m 

relative humidity, percent 

Surface Radiation 

downdraft velocity ratio, unitless 

temperature, °C 
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T s = surface temperature, °C 

UAV = uninhabited air vehicle 

w = updraft velocity, m/s 

w = convective velocity scale, m/s 

Wi = intermediate variable in downdraft velocity calculation, m/s 

w c = updraft velocity corrected for environment sink, m/s 

w D = downdraft velocity, m/s 

w e = environment sink velocity, m/s 

Wpe a k = maximum vertical velocity, m/s 

w s = mixing ratio, g/kg 

W = average updraft velocity, m/s 

X = test area length, m 

Y = test area width, m 

z = aircraft height, m 

z,- = convective mixing-layer thickness, m 

P = Bowen ratio, unitless 

p = moist air density, kg/m 3 

0 O = surface potential temperature, °K 

a = standard deviation 

r 0 = dry adiabatic lapse rate, °C/m 

v|/ = angular position, rad 

0Q = daily average surface potential temperature, °K 


I. Introduction 

O ne relatively unexplored way to improve the range, duration, or cross-country speed of an autonomous aircraft 
is to use buoyant plumes of air found in the lower atmosphere, known as thermals or updrafts. Updrafts occur 
when the air near the ground becomes less dense than the surrounding air because of heating or humidity changes at 
the surface of the Earth. Piloted sailplanes rely solely on this free energy to make flights of more than 
2000 kilometers. 1 Soaring birds such as hawks and vultures have been observed to circle for hours without flapping 
their wings; and Frigatebirds are known to soar continuously, day and night, using updrafts. 2 

Many uninhabited air vehicles (UAVs) have similar sizes, wing loadings, and missions as soaring birds and 
sailplanes. Such missions that could allow small UAVs to take advantage of updrafts encompass remote sensing, 
surveillance, atmospheric research, communications, and Earth science. A previous study using a simple 
UAV simulation with updrafts calculated from measured surface and balloon data found that a 2-hr nominal- 
endurance UAV can potentially gain 12 hr of flight time during the summer and 5 hr of flight time during the winter 
using updrafts. 3 

Researchers have used various updraft models to study optimal soaring paths, 4 glider design, 5 and cloud 
formation. 6 Wharington used a simple updraft model to develop algorithms for a soaring UAV. 7 These models 
generally provide a vertical velocity distribution for a given radius and peak updraft velocity. This paper presents a 
unique updraft model in that statistical values for updraft velocity, convective mixing-layer thickness, updraft 
spacing, and updraft size are determined from measured data. 

II. Convective-Layer Scale Factors 

The convective velocity scale, w*, and convective mixing-layer thickness, z,-, were calculated from surface and 
rawinsonde balloon measurements taken at the National Oceanic and Atmospheric Administration (NOAA) Surface 
Radiation (SURFRAD) station in Desert Rock, Nevada (lat. 36.63 deg N, long. 116.02 deg W, elev. 1007 m). The 
convective layer is the lowest region of the atmosphere where significant mixing occurs. During calm conditions, 
buoyant plumes of air that have been heated at the surface cause local mixing. These plumes of air, called updrafts 
or thermals , can have significant upward velocity. 

SURFRAD instrumentation includes a radiometer platform, meteorology tower, and solar tracker. Measurements 
were taken every 3 min. A rawinsonde station is also collocated with the SURFRAD site where rawinsonde balloons 
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were launched every 12 hr. This study used surface temperature, wind, and radiation measurements as well as 
balloon-measured temperature and humidity collected from the entire year of 2002. 8 

The convective mixing-layer thickness, z,-, as Fig. 1 shows, is the maximum height-above-ground that updrafts 
generally obtain. The mixing-layer thickness was calculated using predawn rawinsonde balloon data and measured 
surface temperatures. The convective mixing- layer thickness was estimated by finding the intersection of the 
predawn balloon-measured temperature profile as a function of altitude with the line given by Eq. (1). 



Figure 1. Simplified representation of atmosphere with updrafts showing convective mixing-layer 
thickness, z ;• 


t = -r D h + t s (i) 

The line defined by Eq. (1) gives a temperature, T, for every height, z, using the dry adiabatic lapse rate, r fl , of 
0.00975 °C/m. Figure 2 shows an example of the z, calculation. 



Temperature, °C 


Figure 2. Example of Zi calculation. 

A surface heat budget was used to calculate the convective velocity scale, w . The first step in this process is to 
calculate the sensible heat flux, Q H , given by Stull in Eq. (2). 9 
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In this equation, Q H is the sensible heat flux, Q s is the net radiation at the surface, Qg is the heat into the ground, 
and P is the Bowen ratio, defined as the ratio of sensible to latent heat fluxes at the surface. This study uses a Bowen 
ratio of 5, given for semi-arid regions. 9 The net radiation, Q s , was measured directly by the SURFRAD station; and 
the heat into the ground, Qg, was calculated applying Eq. (3) taken from Stull 9 using the percentage method for 
daytime calculations. 


Qg = 01 Q s 

Sensible heat flux was then converted to its kinematic form using Eq. (4). 


(3) 


Qh = 


Qh 
P Cp 


(4) 


Values for the density of moist air, p, and the specific heat of dry air, C p , used in Eq. (4), were 1.210 kg/nr 1 and 
1004.67 J/kg°K, respectively. The surface virtual potential temperature flux, Q ov , was then calculated with Eq. (5) 
using the definition of virtual potential temperature for unsaturated air. 9 


Qov =0«( 1 + o.6iwO 


(5) 


Equation (5) was solved using the mixing ratio, w s , calculated from Eq. (6). 

622ev 

= (6) 

p-ev 

In Eq. (6), p is the surface pressure and ev is the vapor pressure given by Eq. (7). 

ev = rh'esl 100 (7) 

Equation (7) was solved using the measured relative humidity, rh, and saturated vapor pressure, es, calculated 
from Eq. (8). 


es = 6.1 12e 


17.672; 
2-, +243.5 


( 8 ) 


Equation (8) was solved using measured surface temperature, T s . Equations (2) to (8) determine the surface heat 
budget required to calculate Qov- Equation (9) calculates the convective scaling velocity. 
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In Eq. (9), g is the gravitational constant and 0 q is the daily average surface potential temperature given by 
Eq. (10). 
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Equation (10) was solved using measured surface temperature, T s and pressure,/!, and a reference pressure, p 0 , 
of 1000 mb. The sample index,/' = 1, refers to sunrise and j = n refers to sunset. The convective velocity scale, w , is 
used primarily to calculate the updraft vertical velocity. The convective scale velocity was set to zero if the surface 
wind velocity was greater than 12.87 m/s (25 knots) to account for the disruptive effect of high winds on updrafts. 

III. Test Point Matrix 

This study uses statistical tools to summarize yearly and daily changes in z, and w into a set of representative 
cases. The histogram in Fig. 3 shows the distribution of w for all daylight hours during the year 2002. Data points 
that produced a convective velocity scale of zero comprise 25 percent of the total set of recorded data and were not 
included in the Fig. 3 histogram. Zero convective mixing-layer thickness, zero net radiation, or high winds can cause 
zero scale velocity. Of these causes, the most prevalent was found to be zero net radiation because of low sun angles 
during early morning and late evening. 
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Figure 3. The w" histogram for all samples during daylight hours of 2002, when updrafts were present. 

A statistical representation of w* was found using percentiles. The following percentiles were chosen as test 
points: 2.3, 15.9, 50.0, 84.1, and 97.7. The five test cases are referred to as -2o, -la, mean, +la, and +2a because 
they are analogous to the standard deviation and mean values of a Gaussian probability density function (PDF). The 
use of percentiles allowed the calculation of a simple set of test cases to represent the distribution of convective 
velocity scale. 

A distribution of all z, values for each selected w in the test point matrix was created by collecting all z, points 
corresponding to w values that fell within a window of ±0.1 m/s from the selected test point w . Figure 4 shows an 
example of the z, data selection for a given w . A Gamma PDF was used to model the statistics of z,-. The Gamma 
PDF coefficients that best fit the data were used to calculate the -la, mean, and ±1 a points in the test point matrix. 
Table 1 gives the resulting set of convective scale factor test points. 
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Figure 4. Convective mixing-layer thickness, zt, as function of convective velocity scale, w . (Dashed lines 
indicate z ,■ window for the example case of w = 4.08 m/s.) 


Table 1. Convective scale test points during times when updrafts are present. 


Description 

w * , m/s 

-la Zj, m 

mean z h m 

±1gz„ m 

-2a w 

0.46 

25.6 

53.6 

97.4 

1 * 
-law 

1.27 

150 

210 

1007 

mean w* 

2.56 

767 

1401 

2319 

±la w 

4.08 

2134 

2819 

3638 

±2 a w 

5.02 

2913 

3647 

4495 


The test conditions given in Table 1 indicate the year-round statistical properties of convective lift during times 
when updrafts were present. These test points include variations in convection as a result of time of day and time of 
year. Table 1 shows large variations of w and z„ illustrating the challenge presented to a soaring aircraft that must 
use as many updraft sizes and strengths as possible. 

Table 2 presents monthly trends. These points represent the mean and maximum values of w during the daylight 
hours of the selected month. The z, test points in Table 2 were found by taking the mean of all z,- points in the 
selected month that fell within a window of ±0.1 m/s from the selected test point w . Test points given in Table 2 
reveal the seasonal variations in updraft strength and altitude. Daily variations are generally sinusoidal, with zero w 
at sunrise, maximum w near noon, and zero w at sunset. 
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Table 2. Monthly convective scale test points for all times between sunrise and sunset. 



Jan 

Feb 

Mar 

Apr 

May 

Jun 

Jul 

Aug 

Sep 

Oct 

Nov 

Dec 

mean w*, m/s 

1.14 

1.48 

1.64 

1.97 

2.53 

2.38 

2.69 

2.44 

2.25 

1.79 

1.31 

1.26 

Zj for mean w*, m 

504 

666 

851 

1213 

1887 

1728 

1975 

1755 

1382 

893 

627 

441 

max w \ m/s 

3.59 

3.97 

4.89 

5.53 

5.49 

5.51 

6.30 

5.64 

5.97 

4.57 

4.55 

4.11 

Zj for max w , m 

1800 

1970 

3900 

2380 

3833 

4027 

3962 

4940 

2460 

3285 

1783 

1680 


IV. Updraft Calculations 

Each w* and z t scaling parameter given in Table 1 and Table 2 can be used to calculate updraft velocity, radius, 
and vertical velocity distribution using the equations provided in this section. Equation (11) taken from Lenschow 10 
calculates average updraft velocity using height, z, and the convective-layer scale parameters, w and z h 
from Table 1 or Table 2. 
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( 11 ) 


Figure 5 is a plot of Eq. (11). The maximum vertical velocity of 0.45 w occurs at a height ratio, z/z,-, of 
approximately 0.25. Also note that the updraft velocity is negative for height ratios greater than 0.9. The height at 
which peak vertical velocity occurs is lower than those published by Bradbury. 11 The dissimilarities may occur 
because of differing climate, topography, or definition of the convective mixing-layer thickness. 
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Figure 5. Updraft mean vertical velocity ratio as function of height ratio. 
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Equation (12) calculates updraft outer radius, r 2 , as determined by Lenschow. 10 
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Equation (12) predicts small updraft outer radius near the ground and an increasingly larger updraft radius as 
height increases. Figure 6 shows updraft outer radius ratio, rVz,-. Updraft radius is independent of velocity in this 
model. 
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Figure 6. Ratio of outer radius to convective mixing-layer thickness, r 2 /Zi, as function of height ratio, zJz.i- 

The distribution of vertical velocity within the updraft varies considerably between updrafts. The literature 
contains various updraft vertical velocity distributions. Wharington 7 ' 12 used Gaussian-shaped updrafts to study 
autonomous soaring techniques and Metzger and Hedrick 4 used simple square-edged lift regions to optimize flight 
paths for cross-country flight. Many references depict updrafts to be quite random in size and updraft velocity 
distribution. In the book Cross Country Soaring , 13 Helmut Reichmann states that thermals are “seldom, if ever, 
round” and “the thermals that depart from the norm are most likely the norm itself.” 

Flight test results were used in an attempt to model the distribution of vertical velocity for the updraft model 
presented in this paper. Flight tests conducted by Konovalov 14 show two basic vertical velocity profiles, called 
type-a and type-b. Type-a updrafts are strong and have a wide area of nearly constant lift. Type-b updrafts are 
weaker and have increasingly more lift toward the center. Konovalov also correlated the frequency of occurrence of 
type-a and type-b updrafts to the updraft diameter. This paper uses Konovalov’s data to correlate the vertical 
velocity distribution of an updraft to its outer radius. Figure 7 shows a revolved trapezoid that was used to fit 
Konovalov’s type-a and type-b updrafts by adjusting the inner radius, r,. This study uses a revolved trapezoid 
because the trapezoid provided a good approximation to the vertical velocity profiles given by Konovalov. A fit of 
Konovalov’s data yields the piecewise function for the updraft radius ratio, r I /r 2 , given in Eq. (13). 
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Figure 7. Revolved trapezoid vertical velocity distribution. 


fO.OO 1 1 * ?2 +0.14 for /"2 < 600m 
[0.8 else 


(13) 


The average vertical velocity of a revolved trapezoid updraft distribution is found by dividing the volume of the 
revolved trapezoid shape by the base area. Equation (14) gives the resulting relationship. 
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Evaluating Eq. (14) and solving for w peak gives the relationship found in Eq. (15). 

3 w(r 2 3 -r 2 2 r 7 ) 
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(15) 


Equation (15) determines the peak value, w pea k, of a revolved trapezoid updraft distribution for a given average 
updraft velocity, w T . The vertical velocity distribution was further refined by fitting a bell shape given in Eq. (16) to 
approximate the revolved trapezoid distribution. A bell shape was chosen because the bell shape defines a smooth 
vertical velocity distribution. 
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In Eq. (16), k^ 4 are shape constants and w D is the downdraft velocity. This equation describes a family of curves 
designed to fit the trapezoid shape given by Eq. (13) for a range of discrete r t lr 2 values. Seven evenly spaced points 
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were chosen to represent the range of possible vertical velocity distributions. Table 3 gives the shape constants used 
in Eq. (16) for each rjlr 2 ratio. Figure 8 shows the resulting family of updraft velocity distributions. The downdraft 
velocity term, w D , simulates the toroid-like downward velocity found on the outer edge of updrafts as described by 
Bradbury. Figure 9 illustrates the toroid shape of the updraft as the updraft increases in height. Eqs. (17) and (18) 
calculate the downdraft velocity. 


Table 3. Shape constants for bell-shaped vertical velocity distribution. 


r l /r 2 

ki 

k 2 

k 2 

k 4 

0.14 

1.5352 

2.5826 

-0.0113 

0.0008 

0.25 

1.5265 

3.6054 

-0.0176 

0.0005 

0.36 

1.4866 

4.8354 

-0.0320 

0.0001 

0.47 

1.2042 

7.7904 

0.0848 

0.0001 

0.58 

0.8816 

13.972 

0.3404 

0.0001 

0.69 

0.7067 

23.994 

0.5689 

0.0002 

0.80 

0.6189 

42.797 

0.7157 

0.0001 



Position, r/r-, 


r l/ r 2 = 0-14 
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Figure 8. Bell-shaped vertical velocity distributions. 
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Figure 9. Development of updraft toroid shape as height increases. 
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for 0.5 < — < 0.9 


else 


(18) 


Additionally, Eqs. (17) and (18) modify the vertical velocity shape given in Eq. (16) to have zero net vertical 
velocity when z/z ; - = 0.9. The height ratio 0.9 was chosen because the vertical velocity of the updraft is negative at 
altitudes above 0.9 as given in Eq. (11). Figure 10 gives vertical velocity profiles for low and high updrafts. The 
vertical velocity profiles show increasing downdraft at the outer edge of the updraft as height is increased. Figure 10 
was generated with the convective velocity scale, w , and convective mixing-layer thickness, z ; , of 2.56 m/s and 
1401 m, respectively. 
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Figure 10. Vertical velocity profiles of a typical updraft with increasing height. 

V. Updraft Spacing and Environment Sink 

Updraft spacing was determined from Lenschow 10 for a constant z/z, of 0.4. Equation (19) gives the resulting 
relationship. 


N l z, 

L 


= 1.2 


(19) 


In Eq. (19), N L is the number of updrafts encountered along a straight line of length, L. Equation (20) solves for 
the number of updrafts within a given area of length, X, and width, Y, and a given updraft outer radius, r 2 . 


N = 


0.6 YX 


~i r 2 


( 20 ) 


In a given test area, all upward moving air in the form of updrafts is balanced by an equal amount of downward 
moving air, known as environment sink, in the surrounding atmosphere. Environment sink is applied whenever the 
aircraft is not inside an updraft. Conservation of mass was used to determine the environment sink velocity for a 
given test area. Equation (21) gives the resulting relationship. 
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Equation (21) relates the environment sink 
downdraft velocity ratio, s w j, given in Eq. (22). 


velocity, w e , to the average updraft velocity, 


W T , using the 


s wcl 
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( 22 ) 
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The downdraft velocity ratio is used to modify the environment sink calculation to account for the downward 
moving air defined in Eqs. (17) and (18). The vertical velocity distribution of the updraft can be blended to match 
the environment sink at the outer radius of the updraft using Eq. (23). 
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1 - 




+ w„ 
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(23) 


peak 
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Equation (23) stretches the vertical velocity distribution to maintain the maximum vertical velocity at the center 
of the updraft while allowing a smooth transition to the environment sink at the outer edge of the updraft. 


VI. Application 

This paper presents an updraft model for a specific intention — to use the model during the design and simulation 
of autonomous soaring UAVs. This section introduces an example set of updraft calculations for a typical simulation 
case. Appendix A gives a MATLAB® (The MathWorks, Natick, Massachusetts) script for this example. The 
convective scale factors, w and z„ were taken from the mean values given in Table 1 to be 2.56 m/s and 1401 m, 
respectively, for this example. The height, z, used in this example was chosen as 280 m, giving a height ratio, z/z,-, of 
0.2. Calculating the updraft outer radius, r 2 , from Eq. (12) yielded 79.4 m. The number of updrafts in the test area 
(defined by X= 1000 m and Y= 1000 m) calculated from Eq. (20) was five, using r 2 = 79.4 m and z t = 1401 m. The 
resulting five updrafts can be positioned anywhere within the test area. Updrafts were positioned evenly along a 
diagonal line in this example for clarity. The vertical velocity of the test area was calculated with Eqs. (11) to (18) 
and (21) to (23) of this paper for every point within the test area using a grid spacing of 10 m. Figure 1 1 shows the 
resulting vertical velocities of the test area. Figure 12 shows the vertical velocity distribution of a single updraft as a 
check case for the code given in Appendices A and B. Figure 12 of this paper was produced with the MATLAB 
script given in Appendix A. Electronic copies of this code can be obtained from the author. 



Figure 11. Example test area (calculated from MATLAB code given in Appendix A). Updrafts are 
positioned along a line for clarity. 
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Y position, m 


Figure 12. Vertical velocity profile check case (calculated from MATLAB code given in Appendix A). 

The updraft positions should be randomly chosen during an actual application. Updraft positions should be held 
for 20 min of simulation time and then re-calculated to account for the duration of a typical updraft. Updraft 
duration can range from 5 to 30 min, indicating that this parameter should be varied during simulation testing of 
autonomous soaring algorithms. 


VII. Future Work 

The updraft model presented in this paper is useful for the development of autonomous soaring guidance and 
control but does not model all of the significant characteristics of naturally occurring updrafts. Naturally occurring 
updrafts have time dependant vertical velocities and radiuses. Additionally, naturally occurring updrafts drift and 
change shape when wind is present; and naturally occurring updrafts merge together to form larger and stronger 
updrafts as height increases. Downdrafts are also known to form along with updrafts during times of strong 
convection. Future updraft models, therefore, should include these effects as well as include convective scale 
parameters for locations and topographies other than the desert. 

VIII. Concluding Remarks 

This paper presented a model of convective updrafts using convective scale parameters calculated from ground 
and balloon measurements taken at Desert Rock, Nevada. Convective scale velocity and convective mixing layer 
thickness test cases were produced using statistical tools. Statistical representation did not include measurement data 
during times of zero convection, representing 25 percent of the total measurement set. A statistical representation of 
the convective mixing-layer thickness values corresponding to the set of convective scale velocities was produced 
by fitting the data to a Gamma probability density function. Convective scale parameters were then used to calculate 
mean updraft vertical velocity and radius using equations from Lenschow. The vertical velocity distribution within 
the updraft was calculated by fitting the data in Konovalov to a family of bell-shaped curves. The toroid structure of 
updrafts as height increases, described in Bradbury, was included in the vertical velocity distribution of this model. 
Updraft spacing was calculated from a relationship given in Lenschow; and environment sink was calculated using 
the conservation of mass. This paper also presented an example use of the equations within the paper using the 
MATLAB code provided in the appendices. Taken together, the equations given herein describe an updraft model 
that is suitable for the preliminary design and simulation of guidance and control for soaring uninhabited air 
vehicles. 
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Appendix A. MATLAB Code for Example Application 

%This script will run a check case of the NASA DFRC updraft model. 

O, 
o 

%Michael J. Allen 

%NASA Dryden Flight Research Center 
%2005 

O, 
o 

%DEFINE UPDRAFT PARAMETERS 
wstar=2 .56; 
zi=1401; 
z=280 ; 
wgain=l ; 
rgain=l ; 

% DEFINE AREA 
X= 1 0 0 0 ; 

Y=1000 ; 

A=X* Y ; 

%CALCULATE OUTER RADIUS 
zzi=z/ zi; 

r2= ( . 102*zzi^ (1/3) ) * (1- ( . 25*zzi) ) *zi; %Equation 12 from updraft paper 

%CALCULATE NUMBER OF UPDRAFTS IN GIVEN AREA 
N=round ( . 6*Y*X/ (zi*r2) ) ; 

%SET PERTURBATION GAINS FOR EACH UPDRAFT 

wgain ( 1 : N) =1 ; %multiplier on vertical velocity 

rgain ( 1 : N) =1 ; %multiplier on radius 

%PLACE UPDRAFTS IN A LINE (NORMALLY THIS IS DONE RANDOMLY) 
for kn=l:N, %for each updraft 

xt (kn) =kn*X/ (N+l) ; 
yt (kn) =kn*Y/ (N+l) ; 

end 

%DEFINE GRID OF TEST LOCATIONS 
xc= [0 : 10 :X] ; 
yc= [0 : 10 : Y] ; 

xx=xc ' *ones ( 1 , length (yc) ) ; 
yy=ones (length (xc) , 1) *yc; 
zz=ones (size (xx) ) . *z; 
for kx=l : length (xc) , 

disp ( [num2str (kx) 'of ' num2str 
for ky=l : length (yc) , 

%CALL UPDRAFT FUNCTION 

w (kx, ky) =run_model2_3 (xx (kx, ky) , yy (kx, ky) , zz (kx, ky) , ... 
xt, yt,wstar, wgain, rgain, zi , A, 1 ) ; 

end 

end 

%PLOT UPDRAFT FIELD. THIS WILL CREATE FIGURE 11 
figure; orient tall 
mesh (xx, yy, w) 


%create mat 
%create mat 
%create mat 

(length (xx) ) ] ) 


rix 

of 

X 

values 

rix 

of 

Y 

values 

rix 

of 

z 

values 


%selected w*, m/s 
%selected zi, m 
%test height, m 

%multiplier on vertical velocity 
%multiplier on radius 

%length of test area, m 
%width of test area, m 
%test area, m^2 
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xlabel('X position, m') 
ylabel('Y position, m') 
zlabel ( ' w, m/s') 
set (gca, 'View' , [-22 64] ) 

%PLOT CROSS -SEC I ON OF 1ST UPDRAFT. THIS WILL CREATE FIGURE 12 

figure; orient tall 

hl=plot (yy ( 18 , 1 : 40 ) , w ( 18 , 1 : 4 0 ) , ' k- ' ) ; grid 

set (gca, ' font size ' , 12) 

set (hi, ' linewidth ' , 2) 

xlabel('y position, m') 

ylabel ( ' w, m/s') 


Appendix B. MATLAB Code for Updraft Model 

function [w,r2,wc] = ... 

run_model2_3 (x, y, z, xt, yt , wstar , wgain, rgain, zi , A, sf lag) 

%function [w,r2,wc] = run_model2_3 (x, y, z , xt , yt , wstar , wgain, rgain, zi , A) 

a. 

o 


%Input : 

X 

= Aircraft x position (m) 

o, 

o 

y 

= Aircraft y position (m) 

o, 

o 

z 

= Aircraft height above ground (m) 

o_ 

o 

xt 

= Vector of updraft x positions (m) 

o_ 

0 

yt 

= Vector of updraft y positions (m) 

o, 

o 

wstar 

= updraft strength scale factor, (m/s) 

o, 

o 

wgain 

= Vector of perturbations from wstar (multiplier) 

o, 

o 

rgain 

= Vector of updraft radius perturbations from average 

o_ 

o 


(multiplier) 

o_ 

o 

zi 

= updraft height (m) 

o, 

o 

A 

= Area of test space 

o, 

o 

sf lag 

= 0=no sink outside of thermals, l=sink 


o_ 

o 


%Output : w = updraft vertical velocity (m/s) 

% r2 = outer updraft radius, m 

% wc = updraft velocity at center of thermal, m/s 

o, 

o 

o, 

o 

%Michael J. Allen, NASA DFRC, 2005 
%DEFINE UPDRAFT SHAPE FACTORS 

rlr2shape = [0.1400 0.2500 0.3600 0.4700 0.5800 0.6900 0.8000] ' ; 
Kshape = ... 


[1 

. 5352 

2 . 

5826 

-0 . 

. 0113 

-0 . 

.1950 

0 . 

. 0008; . 

1 

. 5265 

3 . 

6054 

-0 . 

. 0176 

-0 . 

. 1265 

0 . 

. 0005; . 

1 

.4866 

4 . 

8356 

-0 . 

. 0320 

-0 . 

. 0818 

0 . 

. 0001; . 

1 

.2042 

7 . 

7904 

0 . 

. 0848 

-0 . 

.0445 

0 . 

. 0001; . 

0 

. 8816 

13 . 

9720 

0 . 

. 3404 

-0 . 

. 0216 

0 . 

. 0001; . 

0 

. 7067 

23 . 

9940 

0 . 

.5689 

-0 . 

. 0099 

0 . 

. 0002; . 

0 

. 6189 

42 . 

7965 

0 . 

. 7157 

-0 . 

. 0033 

0 . 

. 0001] ; 


%CALCULATE DISTANCE TO EACH UPDRAFT 
N=length (xt ) ; 
for k=l : N, 

xdsq= (x-xt (k) ) ^2; 
ydsq= (y-yt (k) ) ~2; 
dist (k) =sqrt (xdsq+ydsq) ; 

end 
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%CALCULATE AVERAGE UPDRAFT SIZE 
zzi=z/ zi; 

rbar= ( . 102*zzi A (1/3) ) * (1- ( . 25*zzi) ) *zi; 

%CALCULATE AVERAGE UPDRAFT STRENGTH 
wtbar= ( zzi A (1/3) ) * (1-1 . 1 *zzi ) *wstar ; 

%USE NEAREST UPDRAFT 

upused=f ind (dist==min (dist) ) ; 

if length (upused) >1 ; upused=upused ( 1 ) ; end 

%CALCULATE INNER AND OUTER RADIUS OF ROTATED TRAPEZOID UPDRAFT 
r2=rbar*rgain (upused) ; %multiply by random perturbation gain 
if r2<10, %limit small updrafts to 20m diameter 

r2=10 ; 

end 

if r2<600, 

rlr2=. 0011*r2+ . 14; 
else , 

rlr2= . 8 ; 

end 

rl=rlr2*r2 ; 

%MULTIPLY AVERAGE UPDRAFT STRENGTH BY WGAIN FOR THIS UPDRAFT 
wt=wtbar*wgain (upused) ; %add random perturbation 

%CALCULATE STRENGTH AT CENTER OF ROTATED TRAPEZOID UPDRAFT 
wc= (3*wt* ( (r2 A 3) - (r2 A 2) *rl) ) / ( (r2 A 3) - (rl A 3) ) ; 

%CALCULATE UPDRAFT VELOCITY 
r=dist (upused) ; 
rr2=r/r2; %r/r2 

if z<zi, %if you are below the BL height 

if rlr2< . 5* (rlr2shape (1) +rlr2shape (2) ) , %pick shape 

ka=Kshape (1,1) ; 
kb=Kshape (1,2) ; 
kc=Kshape (1,3) ; 
kd=Kshape (1,4) ; 

elseif rlr2<.5* (rlr2 shape (2 ) +rlr2 shape ( 3 ) ) , 
ka=Kshape (2,1) ; 
kb=Kshape (2,2) ; 
kc=Kshape (2,3) ; 
kd=Kshape (2,4) ; 

elseif rlr2<.5* (rlr2 shape ( 3 ) +rlr2 shape ( 4 ) ) , 
ka=Kshape (3,1) ; 
kb=Kshape (3,2) ; 
kc=Kshape (3,3) ; 
kd=Kshape (3,4) ; 

elseif rlr2<.5* (rlr2 shape ( 4 ) +rlr2 shape ( 5 ) ) , 
ka=Kshape (4,1) ; 
kb=Kshape (4,2) ; 
kc=Kshape (4,3) ; 
kd=Kshape (4,4) ; 

elseif rlr2<.5* (rlr2 shape (5) +rlr2 shape (6) ) , 
ka=Kshape (5,1) ; 
kb=Kshape (5,2) ; 
kc=Kshape (5,3) ; 
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kd=Kshape (5,4) ; 

elseif rlr2<.5* (rlr2 shape ( 6 ) +rlr2 shape ( 7 ) ) , 
ka=Kshape (6,1); 
kb=Kshape (6,2) ; 
kc=Kshape (6,3) ; 
kd=Kshape (6,4) ; 
else , 

ka=Kshape (7,1); 
kb=Kshape (7,2) ; 
kc=Kshape (7,3) ; 
kd=Kshape (7,4) ; 

end 


in=rr2 ; 


%CALCULATE SMOOTH VERTICAL VELOCITY DISTRIBUTION 
ws = ( 1 . / ( 1+ ( ka . *abs ( in+kc) ) . A kb) ) +kd*in; 
ws (find (ws<0) ) =0; %no neg updrafts 

else , 

ws = 0 ; 

end %end if low enough 

%CALCULATE DOWNDRAFT VELOCITY AT THE EDGE OF THE UPDRAFT 


if dist (upused) >rl & rr2<2, 
wl= (pi/6) *sin (pi*rr2) ; 
else , 

wl = 0 ; 

end 

if zzi>.5 & zzi<=.9, 
swd=2 .5* (zzi- . 5) ; 

wd=swd*wl ; 
wd (find (wd>0) ) = 0; 
else , 

swd=0 ; 
wd=0 ; 

end 

w2=ws*wc+wd*wt ; 


%downdraft, positive up 


%scale factor for wd for zzi 
%used again later 


%scale updraft to actual 
%velocity 


%CALCULATE ENVIRONMENT SINK VELOCITY 

At=N*pi* (rbar^2 ) ; %total area taken by updraft 

if At>A, error ('Area of test space is too small'); return; end 

if sflag, 

we=- (At*wtbar* (1-swd) ) / (A-At) ; 

we (find (we>0 ) ) =0 ; 
else , 

we=0 ; 

end 

%STRETCH UPDRAFT TO BLEND WITH SINK AT EDGE 
if dist (upused) >rl ; %if you are outside the core 

w=w2* (1-we/wc) +we; %stretch 

else , 

w=w2 ; 

end 


%environment sink rate, 
%positive up (m/s) 

%don't allow positive sink 
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